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We have developed an efficient Monte Carlo algorithm, which accelerates slow Monte Carlo dy- 
namics in quasi-one-dimensional Ising spin systems. The loop algorithm of the quantum Monte 
Carlo method is applied to the classical spin models with highly anisotropic exchange interactions. 
Both correlation time and real CPU time are reduced drastically. The algorithm is demonstrated in 
the layered triangular-lattice antiferromagnetic Ising model. We have obtained the relation between 
the transition temperature and the exchange interaction parameters, which modifies the result of 
the chain-mean-field theory. 
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The application of the Monte Carlo (MC) method to 
the condensed-matter physics has been successful bridg- 
ing between the experimental study and the theoreti- 
cal study. [H The simulational results are now quantita- 
tively compared with the experimental results. We may 
estimate various physical parameters, predict unknown 
properties, and propose new experiments on real mate- 
rials. However, we encounter a difficulty when we apply 
the MC method to the frustrated systems. The MC dy- 
namics slows down, and it becomes very hard to reach 
the equilibrium states. Since frustration has been recog- 
nized to play an important role in novel effects of many 
materials, we somehow have to overcome this difficulty 
to study new properties, new concepts and new function 
of such materials. 

In this Letter we consider the quasi-one-dimensional 
(Q1D) frustrated spin systems. The magnetic exchange 
interaction of this system is highly anisotropic. The inter- 
action along the c axis is much stronger than those within 
the ab plane: | J c \ ;§> \J a b\- The experimental realizations 
of this model are the ABX3-type compounds. [H, 0, [H, IS 0] 
The lattice structure is the stacked triangular lattice with 
the antiferromagnetic exchange interactions. There are 
two reasons for the slow MC dynamics in this system. 
One is frustration, and the other is the long correlation 
length along the c axis. The single-spin-flip algorithm 
cannot change the states of these correlated clusters. 
Koseki and Matsubara [H, 0, [l(| proposed the cluster- 
heat-bath method, which accelerates the MC dynamics 
in Q1D Ising spin systems. When we update a spin state, 
the transfer matrix is multiplied along the c axis. This 
matrix operation takes a long CPU time. The possi- 
ble size of simulation has been restricted to the system 
with \Jc/Jab\ = 10, 36 x 36 x 360 spins, and 2 x 10 6 MC 
steps. [Ill] Considering that the ratio \ J C / J a b\ in real com- 
pounds is in the order of 100, we need to develop another 
algorithm that improves the simulation efficiency. 

We notice that the similar slow-dynamic situation oc- 
curs in the quantum Monte Carlo (QMC) simulation. 
The d-dimensional quantum system is mapped to the 
[d + l)-dimcnsional classical system, on which the simu- 



lation is performed. The additional dimension is called 
the Trotter direction, and its length is called the Trot- 
ter number. The (d + l)-dimensional classical system 
becomes equivalent to the original d-dimensional quan- 
tum system when the Trotter number is infinite. As the 
Trotter number increases, the correlation length along 
the Trotter direction increases, and the dynamics of the 
simulation slows down. 

The simulation in the Q1D system is equivalent to the 
QMC simulation if we regard the Trotter direction as the 
c axis in the Q1D system. For example, the cluster-heat- 
bath algorithm in the Q1D system is equivalent to the 
transfer-matrix MC method^, [3] in QMC. This is the 
main idea of this paper. We know that the continuous 
imaginary-time loop flip algorithm of QMC [IB, 1|| 17, H} 
is very efficient. Therefore, we apply this QMC algorithm 
to the Q1D simulation. The correlated cluster along the 
c axis is flipped by one update trial. We do not suffer 
from the MC slowing-down due to the long correlation 
length. The algorithm was successfully applied to the 
theoretical analysis on the magneto-electric transitions in 
RbCoBr 3 .[l3] The numerical results quantitatively agree 
with the experimental results. The estimates of the in- 
teraction parameters and proposals of new experiments 
were made possible. 

We consider the transverse-field Ising model in two di- 
mension. The Hamiltonian is written as 
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where a x and a z denote the Pauli spin operators, J de- 
notes the exchange interaction parameter, and T(> 0) 
denotes the transverse field. The indices, j and k, de- 
note the spin location on the two-dimensional real-space 
plane throughout in this paper. The bracket (• • • } de- 
notes the interacting spin pairs. We apply the Suzuki- 
Trotter decomposition[12| and map the quantum system 
Hq to the effective classical system H. c , which is written 
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Here, m denotes the Trotter number, (3 denotes the in- 
verse temperature, and <Xjj = ±1 denotes the Ising spin. 
The index, i, denotes the location along the Trotter di- 
rection throughout in this paper. The first term of this 
effective classical system denotes the exchange interac- 
tion between spins on the same Trotter slice. The second 
term is the exchange interaction between spins at the 
same real-space site with the different (nearest-neighbor) 
Trotter slice. 

The effective classical system can be regarded as the 
QfD spin system: 

Hqid = 2^ I _ '^ ab a i,i a i,k ~ Je / J cr i,j (7 »+l,j 
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Here, i denotes the location along the c axis, j and k 
denote the location on the ab plane, and L c denotes 
the linear size along the c axis. The simulation in the 
Q1D system can be substituted for the QMC simula- 
tion with J, r, and m defined above. The sign of J c is 
positive (ferromagnetic) in this expression. In the case 
when it is negative (antiferromagnetic), we transform it 
to the ferromagnetic one by changing the spin notation 
as o;j — > (— 1 )Vij. 

Let us consider the cluster algorithm of the Q1D sys- 
tem. It is the interpretation of the QMC cluster algo- 
rithm [l8j|, where the cluster is only defined along the 
Trotter direction. We define a cluster using only the J c 
part of the QfD Hamiltonian, and consider the J a b part 
as the molecular field to the cluster. We may regard this 
algorithm as the Swendsen-Wang algorithm [201 ] in one 
dimension [the J c part in Eq. (|5J)] under the molecular 
field [the J a b part in Eq. (jHJ)]. The ergodicity and the 
detailed-balance condition are guaranteed. 

The updating procedure is as follows. First, we se- 
lect one location j on the ab plane, and consider the 
spins along the c axis. We define clusters by connecting 
the neighboring spins (<Jij and c^+i^) with the following 
probability p c : 



Pc = 1 - exp[-2/3J c ] (aij = o-i+i.j), (7) 
Pc = {a. l:j ^ <Ti+i,j). (8) 



Let us number the cluster by /. Second, we calculate the 
molecular field hj for each updating cluster I as 
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where i £ I denotes that i belongs to the cluster J, and 
(k) denotes that cjj^ is interacting with (J%j. Finally, we 
flip the cluster state with the following probability pr. 

Pl ~ exp[2^~7i7] + I ' ^ 

We independently try this flip for each cluster. 

The MC correlation time is reduced by this flip but 
the real CPU time rather increases because we have to 
do the connecting procedures for all spins along the c 
axis. We solve this problem by applying the continuous 
imaginary-time cluster flip algor ithm of QMC. [3 E3, El 
We neglect the discreteness of the spin location along the 
c axis. This approximation is possible when both cluster 
length and L c are very long. 

In the continuous version we focus on the locations of 
the cluster edges. The probability that the spin pair of 
0i 



>i,3 — fj+ij is n °t connected is 

exp[-2/3J c ] =C\ 



where £ c is regarded as the correlation length along the 
c axis. The average cluster size coincide with the cor- 
relation length £ c . If we set L c = L a fc£c, the system 
roughly consists of L\ h correlated clusters. It is known 
that the cluster length obeys the Poisson distribution. 
We generate the Poisson random numbers with the mean 
exp[2/3J c ] and regard them the cluster length. Then, 
we place the cluster edges to the c axis from bottom to 
top. Combining these new-generated cluster edges and 
the already-existing ones, we apply the cluster flip with 
the probability Pj. The procedure is shown in Fig. [I] 

The present continuous c axis version benefits from 
the memory reduce and the CPU time reduce. We do 
not need to memorize all the spin state. Only the loca- 
tions of the cluster edges and the spin state at the bottom 
edge are necessary. The total memory use and the real 
CPU time are proportional to L^ b . Those for the single- 
spin-flip algorithm are proportional to the total number 
of spins, L^ b tlc- The efficiency gain, £ c , becomes expo- 
nentially large at low temperatures. 

We apply the continuous c axis cluster flip algo- 
rithm to the stacked-triangular lattice antiferromagnetic 
Ising model. It is a model system for the ABX3-type 
compounds. 0, 0, 0, @, 0] The Hamiltonian is written as 
follows. 
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FIG. 1: The updating procedure of the continuous c axis 
version. Black (white) rectangles depict the up-state (down- 
state) spin clusters. The new- generated cluster edges are de- 
picted by broken horizontal lines. Brackets depict the clusters 
to be updated. We update each cluster state independently 
using the probability pi. 



where Sij = \<Ji.j is the spin-1/2 Ising spins, and J\(J2) 
denotes the nearest-neighbor (next-nearest-neighbor) ex- 
change interactions within the ab plane. We consider 
the case where both J c and J\ are antifcrromagnetic 
(J c , Ji < 0), and J2 is ferromagnetic (J2 > 0). 

It is known through the theoretical analyses 13, 21. 22. 
23[ that successive magnetic phase transitions occur. The 
low-temperature magnetic structure is the ferrimagnetic 
state. There exists a partially-disordered (PD) phase 
between the paramagnetic phase and the ferrimagnetic 
phase. In the PD phase, one of three sublattices is com- 
pletely disordered, while the other two sublattices take 
antifcrromagnetic configurations. It is considered that 
the phase transition between the paramagnetic phase and 
the PD phase is the second-order transition. We refer to 
the transition temperature as Tni- 

We compare the equilibration and the real CPU time of 
the present algorithm with the results of the single-spin- 
flip algorithm. We set J c = -97.4K, J x = -2.44K, J 2 = 
0.142K, and perform the simulation at T = 25K. The 
system is in the ferrimagnetic phase at this temperature. 
The linear lattice size of the ab plane is set as L a b = 
95. The correlation length along the c axis is roughly 
estimated as £ c ~ cxp[/3|J c |] = 49, and the linear size 
along the c axis is set as L c = L a b£ c = 4655. The effective 
number of spins is more than 42 millions. We observe 
the relaxation functions of the structure factors defined 
as follows. 



\?7=q,/3,7 

fl = ((m a + mp + m 7 ) 2 ) , 



(13) 
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where m a , mp, and m 7 are three sublattice magnetiza- 
tions in the triangular lattice. The 1/3-structure factor, 
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FIG. 2: (Color online) The relaxation function of the struc- 
ture factors, /1/3 and 9/i . 



MC steps [MCS] 


Cluster flip [s] 


Single-spin flip [s] 


100 


39 


568 


200 


77 


1155 


500 


195 


2910 


10000 


3919 


58560 



TABLE I: The real CPU time for each Monte Carlo step is 
compared. The simulations were performed on the Core 2 
Duo E6600 processor at 2.4GHz using the Intel compiler. 



/jy 3 takes a finite value when the ferrimagnetic state or 
the PD state is realized. It detects the phase transition 
between the PD phase and the paramagnetic phase. The 
phase transition between the PD phase and the ferrimag- 
netic phase is detected by the structure factor, fx. 

Figure [2] shows the relaxation functions of both struc- 
ture factors. We start the simulation from the perfect fer- 
rimagnetic state, where structure factors take fy 3 = 1 
and ff = 1/9. The data of two algorithms converge 
to the same value. It guarantees the equilibration of 
the simulation. The cluster algorithm realizes the equi- 
librium state roughly 300 times earlier than the singlc- 
spin-flip algorithm. Table [I] compares the real CPU 
time. The present cluster algorithm achieves the 15 times 
faster simulation. This difference comes from the ratio 
L c / L a b = £ c . 

We focus on the transition temperature between the 
paramagnetic phase and the PD phase, Tni- The chain- 
mean- field theory [2l| gives the relation among Tni, Ji, 
and J c , which is written as 
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Using the present cluster algorithm we estimate Tni for 
various choices of Ji/J c ranging from 0.001 to 0.5 and 
•h/Jc ranging from -0.05 to -0.0015. The behavior of 
Tni with respect to J c , J±, and Ji is obtained. Here, the 
nonequilibrium relaxation method[24j is applied. We ob- 
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FIG. 3: (Color online) Relation between the exchange in- 
teraction parameters and the transition temperature Tni ob- 
tained by the Monte Carlo simulation. Arrows depict the 
point where IJ2I = |Ji|/2 for each choice of JilJ c . The 
numerical results fall onto the single function as long as 
\M < |Ji|/2. 



tain the transition temperature by the behavior of the 
relaxation functions of the structure factor, fy 3 . It con- 
verges to the finite value when the temperature is below 
Tni and decays exponentially when the temperature is 
above Tni- The algebraic decay is exhibited at Tni- 

We find that most of our numerical results are well- 
fitted by the following expression. 
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Only the coefficient of Ji differs from the chain-mean- 
field result. The change of the coefficient can be regarded 
as the reduction of the effective coordination number. j2f| 
The fitting is plotted in Fig. [3] Arrows in the figure 
depicts the data when | = \J\ |/2 for each choice of J 2- 
The data deviate from the relation, Eq. (fT6|) . when | J2I > 
Ji|/2. Since the chain-mean-field relation, Eq. (JT5J) , has 
been used to estimate the interaction parameter from 
the experimental results, the present relation, Eq. (|16[) . 
improves the estimate. 

We have introduced the cluster flip algorithm suitable 
for the quasi-onc-dimensional frustrated Ising spin sys- 
tems. The numerical efficiency is improved as we lower 
the temperature and/or as we increase the anisotropy ra- 
tio, I J c /Jab\- Other algorithms mostly fail in this situa- 
tion. The realistic simulations (or emulations) for real 
compounds are made possible. The quantitative MC 
analyses to the experimental results may help develop- 
ments in the material science. Simulations under the 
magnetic field is possible. We may include the field term 



into the molecular field term, hi. 
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